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We study the folding of the regular two-dimensional triangular lattice embedded in the 
regular three-dimensional Face-Centred Cubic lattice, a discrete model for the crumpling 
of membranes. Possible folds are complete planar folds, folds with the angle of a regular 
tetrahedron (71°) or with that of a regular octahedron (109°). We study this model in 
the presence of a negative bending rigidity K, which favours the folding process. We use 
both a cluster variation method (CVM) approximation and a transfer matrix approach. 
The system is shown to undergo two separate geometrical transitions with increasing \K\: 
a first discontinuous transition separates a phase where the triangular lattice is preferen- 
tially wrapped around octahedra from a phase where it is preferentially wrapped around 
tetrahedra. A second continuous transition separates this latter phase from a phase of 
complete folding of the lattice on top of a single triangle. 
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1. Introduction 



The statistical mechanics of membranes and flexible two-dimensional interfaces is an 
extremely rich and rather challenging field Of fundamental importance is the phase 

diagram of such systems and the precise nature of the remarkable orientationally ordered 
phases that they possess. Faced with a difficult problem it is often useful to explore 
extreme limits. One such limit is that of a flexible phantom fixed-connectivity membrane 
with infinite elastic stiffness. In this limit the in-plane "phonon" degrees of freedom are 
frozen, leaving only out-of-plane "bending" modes. In a 1/d expansion (where d is the 
dimension of the embedding space) such a membrane exhibits a continuous crumpling 
transition at a positive critical bending rigidity . 

Discretising this type of membrane as a triangular lattice leads to a model of the 
statistical mechanics of folding. Finally one can discretise the space of allowed local folds 
to obtain a completely discrete "vertex" model. The simplest case is planar folding, in 
which the local folds are no- folds or complete folds This corresponds to folding the 

triangular lattice onto itself. Planar folding admits many equivalent representations, each 
of them interesting. It may be formulated as a 11-vertex model equivalent to the three- 
colouring problem of a triangular lattice [|5]||. For zero bending rigidity this model was 
solved by Baxter, in its dual formulation, as the combinatorics of edge-colourings of the 
honeycomb lattice [jTj. It may also be mapped to the fully packed loop model with fugacity 
2 describing non-intersecting loops on a honeycomb lattice such that every vertex is covered 
by a loop |fj§. Finally an interpretation of the Z3 colour variable as a Potts spin living 
on the centres of the bonds maps planar folding to a description of the ground state of the 
three-state antiferromagnetic Potts model on the Kagome lattice |jiO| |. Adding a positive 



bending rigidity to the planar folding model one finds a first-order phase transition to a 
completely unfolded phase [pUTTf. 

Since planar folding is very restrictive, it is natural to consider enlarging the discrete 
space of allowed folds. In reference [|12| planar folding was generalised to a particular three- 
dimensional folding allowing four types of discrete folds. This corresponds to choosing 
the embedding space to be a three-dimensional Face-Centred Cubic lattice. This 3-d 
folding model has a formulation as a 96-vertex statistical mechanics model and, once 
again, undergoes a first-order transition to an ordered phase at a definite positive bending 
rigidity |13] , |14f . It also admits a formulation as a dressed three-colouring problem |L2[ . 



In this paper we extend the study of the 96-vertex model to negative bending rigidity. 
In such a model folding is energetically favoured. Negative bending rigidity might in 
principle be realisable in a polymerised membrane with strong next-to-nearest neighbour 
interactions [15] or for a membrane in a sufficiently bad solvent. 

Utilising both the cluster variation method and a transfer matrix formalism we find 
two distinct geometrical phase transitions at negative bending rigidity. As K becomes more 
negative starting from zero, there is a first-order transition from a phase of preferential 
octahedral wrapping to a phase with preferential tetrahedral wrapping, followed by a 
continuous piling transition to a single plaquette. We determine a set of three order 
parameters that completely characterise the separate phases observed. 

The outline of the paper is as follows. In Section 2 we review the construction of the 
96-vertex model. In Section 3 we discuss four natural submodels of the 96-vertex model — 
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these correspond to different restrictions on the allowed folds. We also introduce a duality 
transformation which maps any allowed folding configuration to another allowed folding 
configuration. In Section 4 we study the statistics of folding for negative bending rigidity 
and introduce the three geometrical order parameters. The phase diagram is determined 
by the Cluster Variation Method (CVM) and by a transfer matrix method. In Section 5 
we discuss the incorporation of external fields. Finally we summarise our results and their 
significance in Section 6. 

2. The 96-vertex model 

2.1. General definition of the folding model 

In this section, we first recall the model of discrete folding of the regular 2-dimensional 
triangular lattice embedded in the regular 3-dimensional Face-Centred Cubic (FCC) lattice. 
A detailed derivation of the model can be found in ||12|| . As explained below, the model 
is easily expressible in terms of two Z 2 variables a and z, defined on the faces of the 
triangular lattice. The associated domain walls indicate the different types of folds. The 
spin variables a and z are subject to two basic folding rules, leading to a re-expression of 
the model as a 96-vertex model on the triangular lattice. 




Fig.l: The Face-Centred Cubic lattice viewed as a packing of 3-d space with 
octahedra and tetrahedra. 

A folding of the triangular lattice in the FCC lattice is simply a mapping which assigns 
to each vertex of the triangular lattice a vertex in the FCC lattice, with the constraint 
that neighbouring vertices on the triangular lattice are sent to nearest neighbours in the 
FCC lattice. The FCC lattice can be viewed as consisting of a filling of space by regular 
octahedra complemented by regular tetrahedra (see Fig. 1). Nearest-neighbour triplets 
in the FCC lattice thus form equilateral triangles. In the folding process each elementary 
triangle of the triangular lattice is mapped onto one of these equilateral triangles of the 
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Fig. 2: An example of folding in the FCC lattice of an elementary hexagon 
cut out from the triangular lattice. Nearest neighbours in the hexagon remain 
nearest neighbours after embedding in the FCC lattice. 



acute 




Fig. 3: The four possible folding angles between two adjacent triangles. The 
neighbour of the dark triangle may (i) occupy the symmetric position in the 
same plane (no fold), (ii) be itself on top of the dark triangle (complete fold), 
(iii) lie on the same tetrahedron (i.e. form an acute angle), (iv) lie on the 
same octahedron (i.e. form an obtuse angle). 



FCC lattice. Fig. 2 shows an example of folding of a piece of triangular lattice consisting 
of an elementary hexagon excised from the lattice. 

After folding, two adjacent triangles form some relative angle, i.e. their common link 
serves as a hinge and may be (partially) folded or not. There are four possible relative 
positions for two given neighbouring triangles, as depicted in Fig. 3, measured simply by 
the angle 9 of the corresponding fold: 

(i) 9 = 180° — no fold: the triangles are side by side; 

(ii) 9 = 0° — complete fold: the triangles are on top of each other; 

(iii) 9 = arccos(l/3) ~ 71° — fold with acute angle: the two triangles lie on the same 
tetrahedron, and 

(iv) 9 = arccos(— 1/3) ~ 109° — fold with obtuse angle: the triangles lie on the same 
octahedron. 

It has been shown in [12] that these four possible types of folds are actually organised 
as the domain walls for two Z2 variables a = ±1 and z = ±1 defined on each face of the 
triangular lattice. The relative value of both Z2/Z1 and (J2/G1 for two neighbouring triangles 
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indicate which type of fold they form, with the correspondence displayed in Table I below. 
The domain walls for the z variable are the location of the folds which are either acute 
or obtuse, whereas those for the a variable are the location of the folds which are either 
complete or obtuse. The superposition of these two types of domain wall fixes the folding 
state of all the links, specifying the folding state of the lattice (up to a global orientation). 



Z2/Z1 


0-2/0-1 


angle 


1 


1 


no fold 


1 


-1 


complete fold 


-1 


1 


acute fold 


-1 


-1 


obtuse fold 



Table I: The relative folding state of two neighbouring triangles according 
to their relative values of Zi = ±1 and o~i = ±1 (i = 1, 2). 

Note, of course, that the two global symmetries a — > —a and z — > —z, independently, 
leave the folding invariant. 

2.2. The two basic folding rules 




Fig. 4: The six zi and o~i variables around a given vertex, and the colours q 
of the interior links. 

With the above definitions, the a and z variables on neighbouring triangles are not to- 
tally independent. This can be seen by considering the folding of any elementary hexagon 
made of six neighbouring triangles in the lattice. The relative positions of the six sur- 
rounding vertices of the hexagon can be traced back from the folds of its six internal links. 
After a complete turn around the hexagon, the same absolute position must be found for 
any of these six vertices. This leads to two basic folding rules for the values Oi and Zi 
(i = 1, • • • , 6) on these six triangles (see Fig. 4). 

The first rule involves the variable a only and states that: 

6 

^(T^Omod 3. (2.1) 

i=i 
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This rule is easily understood in terms of 3-colourability of the links of the triangular 
lattice as follows. When embedded in the FCC lattice, the links of the triangular lattice 
must belong to one of the three crystalline planes (200), (020) or (002) which cross the 
octahedra by their three square bases. If we assign a colour to each of these three planar di- 
rections and colour the links of the FCC lattice correspondingly, this induces a 3-colouring 
of the triangular lattice links with three colours (c = 0, 1, 2 mod 3), which, moreover, have 
to be distinct on the three edges of the elementary triangles. One can easily check that 
the colours c% of the six interior links of a hexagon (see Fig. 4) are directly linked to the 
Oi variables by: 



^2 °k mod 3. (2.2) 



k=i+l 



For instance, if <j%I<j\ = — 1, then C2 — cq = a± + o<i = since a complete or an obtuse 
fold puts the link 2 in the same crystalline plane as the link 0. The first basic folding 
rule simply states that the same colour must be found after a complete turn around the 
hexagon, i.e. that cq = c%. 

This first rule on the a variable is not sufficient to ensure that the (zi, Gi) do correspond 
to an actual folded state of the hexagon. This is because the colour only partially specifies 
the actual orientation of the links in 3-dimensional space. A second basic folding rule 
involves both the z and a variables and can be stated as follows: 

{ [ ZiZ i+1 = 1 for c = 0, 1, 2 

iei(c) 

(2-3) 

1(c) = {i : ^^crfc = c mod 3}. 

k=i 

This rule states that the number of changes of the z variable (which indicate folds which 
are either obtuse or acute) occurring on all interior links of a given colour in the hexagon 
must be even. The interpretation of this second rule will be clarified in Section 3 by use of 
a duality transformation. A direct interpretation based on the tetrahedral group can be 
found in |12 . 



It can be shown that the two rules ( |2.1| ) and (|2.3|), applied independently to all the 
hexagons of the original triangular lattice, are necessary and sufficient to define a folded 
state in the FCC lattice. Of course, our definition of folding does not impose any restric- 
tions on possible self-intersections of the triangular lattice, which is therefore considered as 
a phantom object. In the above formulation, the folding problem can thus be understood 
as a constrained system of two Zi variables z and a. 

It is interesting to work out all possible hexagonal configurations around a vertex 
allowed by the two rules ( |2.1| ) and (|2.3| ). If one forgets the global fourfold symmetry 
under global reversal of all a or z and represents only the corresponding domain walls, 
i.e. the four types of folds (i)-(iv) discussed above, one is led to 96 possible hexagonal 
configurations for any of the vertices of the triangular lattice. These are represented in 
Fig. 5. Our folding problem can thus be thought of as a 96-vertex model defined on the 
triangular lattice. 
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Fig. 5: The 96 possible hexagonal configurations around a vertex of the tri- 
angular lattice. Complete folds are represented by thick solid lines, obtuse 
folds by thin solid lines, and acute folds by dashed lines. We have represented 
only the vertices which are not equivalent to each other by a simple rotation. 
We indicate for each vertex the corresponding degeneracy under rotation. 

2.3. The Folding Energy 

For each type of fold, we introduce the folding energy: 

E = K cos(0), (2.4) 

where 6 is the angle of the fold (180° for no fold, 0° for a complete fold, 71° for an acute 
fold and 109° for an obtuse fold). The parameter K is the bending rigidity. In the following 
we shall be interested mainly in the case of negative K, which corresponds to favouring 
folds. In terms of the variables o\, a 2 and z\, z 2 of the two triangles forming the fold, the 
energy simply reads: 

E = -ja 1 a 2 (l + 2z 1 z 2 ). (2.5) 
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The total folding energy is the sum of all elementary folding energies for all links of the 
triangular lattice. 

3. Submodels and duality 

In this section, we first introduce four submodels of the 96-vertex model describing 
four particular restricted folding problems. These models can be seen as particular limiting 
cases of the full 3-d folding problem, where some particular complete ordering is realised 
in the folding process, which correspond to a restriction of the target space to subsets 
of the FCC lattice. We shall see in the next section that similar partial orderings occur 
spontaneously in the unrestricted full 3-d folding problem. We also derive here a duality 
relation for the 96-vertex model. 

3.1. Planar folding: an 11-vertex model 

A first submodel can be obtained by restricting the target space to a single plane 
made of triangles in the FCC, say the crystalline plane (111) which cuts the FCC lattice 
along a regular 2-dimensional triangular lattice. In this case, we are dealing with the pure 
2-dimensional folding problem studied in |4|-|6|]. Here, only complete folds or no-folds are 
allowed. The model can thus be characterised by: 

z = +l (3.1) 

globally on all triangles of the lattice (of course the choice + is arbitrary and a solution 
with z = — 1 globally would be equally acceptable). With this constraint, the second 
folding rule ( |2.3| ) is automatically satisfied and the problem translates into a model with 
a single variable a constrained by the folding rule ( |2.1| ). This model has been studied 
extensively in ||. In particular, the number of folded states of a lattice of N triangles is 
found to grow like q N with q = V3T(l/3) 3 / 2 /2tv ~ 1.208. In terms of vertices, the model 
leads to 11 allowed vertices, which are those of the first row in Fig. 5. 

3.2. Folding on a single octahedron: a 16-vertex model 

One can restrict the target space to a single octahedron in the FCC lattice by simply 
allowing complete and obtuse folds only. Such folds are characterised by a change in 
the value of a on neighbouring triangles. It is thus convenient to introduce a staggered 
variable by dividing the original triangular lattice in the two subsets A and B made of all 
triangles pointing up and all triangles pointing down respectively in the original (unfolded) 
triangular lattice. We define 

cr st = cr on subset A , 
a st = — & on subset B . 
The submodel can then be simply characterised by: 

a s t = +1 (3.3) 
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globally on the triangular lattice. With this constraint, the first folding rule is auto- 
matically satisfied since Y2^=i o"i = on all hexagons. The second folding rule takes a 
particularly simple form since Eq. (|3.3| ) imposes that the interior links of any particular 
hexagon take only two colours, say and 1, which alternate such that 1(0) = {1, 3, 5} and 
1(1) = {2, 4, 6}. The rule ( |2.3|) then translates into a single equation: 

6 

Hzi = l. (3.4) 

i=l 




Fig. 6: The 16 vertices made of complete and obtuse folds only, characterising 
the folding on a single octahedron. 

The two rules (|3.3|) and (|3.4j ) leave us with the 16 vertices of Fig. 6, where the thin 
lines (obtuse folds) indicate the domain walls of the z variable. One can easily show that 
the number of allowed folding configurations for iV triangles grows like q N , with q = \/2. 
Indeed, all the entropy comes from the z degree of freedom and the constraint ( |3.4j ) is 
solved by introducing Z2 variables r\ = ±1 on the vertices of the triangular lattice and 
requiring that the variable z on a given triangle is the product of the three 77 variables 
on the three surrounding vertices. In this formulation the variables r\ are not constrained, 
leading to 2 N ' 2 configurations for the N/2 vertices of the lattice. 

3.3. Folding on a single tetrahedron: an 11-vertex model 

A third submodel corresponds to the folding on a single tetrahedron, meaning only 
complete and acute folds are allowed. The model can be characterised simply by requiring 
that 

z = +a st (3.5) 

for all triangles in the lattice, with a variables satisfying the first folding rule (|2.1|) . Note 
that here again, ( |3.5| ) automatically implies the second folding rule fl2.3|) so that no extra 
constraint has to be imposed on the a variables. In other words, for a variables satisfying 
the first folding rule (|2.1| ), the condition 

\\ (<7Wt)j (t7st)i+l = ] [ (-o-i(Ji+i) = 1 (3.6) 
iei(c) iei(c) 
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is automatically satisfied. This property can be checked directly by inspection of the 22 
sets of <Tj satisfying the first folding rule. In our context, it can be understood simply 
as follows: given an arbitrary hexagonal configuration of a variables around a vertex, 
satisfying the first basic folding rule, this configuration defines a particular planar folded 
state in, say, the (111) plane.0 Given an arbitrary colour c, the image in the FCC lattice 
of the central vertex of the hexagon belongs to exactly one plane of the chosen colour. The 
(111) plane intersects this plane through a line of the chosen colour. Successive triangles 
around the hexagon can lie on the same side of this line or not. The passage from one side 
to the other takes place on a link of the chosen colour c and requires no fold, characterised 
by no change in a, i.e. by a change in cr st . The number of changes in <r s t occurring on 
links of the colour c thus indicates the total number of crossings of the line when making 
a turn around the hexagon. This number must clearly be even since one must get back to 
the same side after one turn. This leads to (|3.6|) . 




Fig. 7: The 11 vertices made of complete and acute folds only, characterising 
the folding on a single tetrahedron. 

The submodel of folding on a single tetrahedron is restricted to the 11 vertices dis- 
played in Fig. 7, made of acute and complete folds only. These vertices are in direct 
correspondence with those of the planar 11-vertex model by replacing no-folds by acute 
folds. This is a particular example of the duality transformation discussed below. 

3. 4- Folding on a single triangle: a 1-vertex model 

Finally, the three submodels above share a common vertex with all folds being com- 
plete, characterised by z = a st = +1 everywhere. This trivial case describes a state of 
complete folding of the triangular lattice on top of a single triangle. 

The above results are summarised in Table II, which lists for each submodel the 
corresponding rules for the a and z variables, the restricted target space and the associated 
number of vertices. 

1 This planar folded state and the tetrahedral folded state will be dual to each other under the 
duality transformation described in the next section. 
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a variable 


z variable 


target space 


vertices 


dual model 


6 

ai = mod 3 

i=l 


z = +l 


single plane 


11 


single tetrahedron 


0"st = +1 


6 
i=l 


single octahedron 


16 


self- dual 


6 

&i = mod 3 

i=l 

0"st = +1 


z = +a st 
z = +1 


single tetrahedron 
single triangle 


11 
1 


single plane 
self- dual 



Table II: Different submodels of the 96-vertex model. 

3.5. The duality transformation 

Given any folded state characterised by its a and z variables satisfying the two basic 
folding rules around all hexagons, we have seen that the first rule automatically implies 
Eq. (|3.6| ). Multiplying this equation with the second folding rule on the z variables, we 
get: 

Y[ {zo- st )i(za st ) i+ i = 1, (3.7) 
iei(c) 

indicating that the global transformation 

(a, z) -»• (a, w = za st ) (3.8) 

maps any folded state of the triangular lattice to another equally acceptable folded state. 
We shall refer to the above transformation as the duality transformation of the 96-vertex 
model. It is clearly an involution, since a^ t = 1, and corresponds geometrically to a global 
interchange of all acute folds and no-folds. 

It is interesting to look at the behavior of the different submodels of the previous 
subsection under the duality transformation. The 16- vertex model is self-dual in the strong 
sense that any of its vertices is invariant under the duality transformation. The 11-vertex 
model describing the planar folding problem and that describing tetrahedral folding are 
dual to each other. In particular the study of reference || on the 2-d folding problem 
applies to the tetrahedral folding problem as well. The duality transformations of the 
different submodels are shown in Table II. 

Finally, let us discuss again the meaning of the second folding rule ( EO)) . By duality, 
this second rule is totally equivalent to (|3.7|). Now ( |3.7| ) simply expresses that the number 
of changes in the variable za st occurring on the links of a given colour must be even. Such 
changes account for obtuse folds or no-folds, which are precisely the two ways of crossing 
the crystalline plane of the given colour. The number of such crossings on a loop must 
clearly be even, implying Eq. (|3.7|) . 

To conclude this section, we would like to mention that the duality transformation is a 
particular case of a more general group transformation for the allowed z configurations at 
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a fixed a global configuration satisfying the first folding rule (|2.1|). Indeed, given two folded 
states (cr, z) and (cr, z') defined on the entire lattice with the same cr, the configuration 
(cr, zz') also defines an allowed folded state, as is clear from the multiplicative nature of 
the second folding rule (|2.3j ). The duality transformation simply corresponds to the allowed 
configuration z' = cr s t- 



4. Geometrical folding transitions 

In this section, we will study the statistics of folding in the case of negative bending 
rigidity K. We will be interested primarily in the three following average values: 

O = (a st ) , 

T = (za st ) , (4.1) 
P=(z), 

with of course — 1 < P, O, T < 1. The average is taken over both all folded configurations 
and all triangles in the lattice. Following Table II, we will refer to P, O and T as the planar, 
octahedral and tetrahedral order parameter respectively, since P = 1 (respectively O = 1 
and T = 1) indicate a situation of pure planar folding (respectively folding on a single 
octahedron and folding on a single tetrahedron). A non-zero value of P will thus be the 
signature that the triangular lattice preferentially stays within planes of the FCC lattice. 
Similarly, we will interpret a non-zero value of O (respectively T) as a tendency for the 
triangular lattice to be preferentially wrapped around octahedra (respectively tetrahedra) . 
We will also consider two-point functions such as: 

02 = (MiiaM , 

T2 = ((^ st ) 1 (z C T st ) 2 ), (4.2) 
P2 = (z lZ2 ) , 

for the correlations on two neighbouring triangles in the lattice. 

4-1. The CVM approximation 

The Cluster Variation Method (CVM) is an approximation of the exact variational 
formulation of the problem |T6| - |T8[1 . Consider our system of variables z, a on the lattice 



and the energy ( |2.5|) summed over all pairs of neighbouring triangles. In terms of a density 
matrix p({cx, z}) for each allowed global folding configuration {cr, z}, the exact free energy 
F of the system can be obtained from: 

F=[j2p(W^})[E({a,z})+\np({a,z})] .. (4.3) 

{cr,z} 
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Here we absorb the temperature factors ksT in the definition of E and F. The subscript 
rain means that the expression is taken at its minimum with respect to p({o~, z}), with the 
normalisation constraint: 

5>({<t,z}) = 1. (4-4) 

{a,z} 

This is the well-known variational principle. Average values are calculated according to: 

where p m i n ({a, z}) is the density at the minimum of the free energy functional (^.3| ). 

The CVM approximation is obtained by truncating the generalised cumulant expan- 
sion for the entropy part of (|4.3| ), S = — J2{a z} P({ a i z }) l n /°({0, z }), to a se t of clusters 
with a fixed maximal size. Here we choose this set to be made of all elementary hexagons 
in the system and all their subclusters, i.e. we take the elementary hexagons, made of six 
neighbouring triangles, as the largest clusters. We thus use the truncated entropy: 

scvm= J2 Sh ~ E Sp+ E St - ( 4 - 6 ) 

hexagons H pairs P triangles T 

Here we call Sh = ~Ph ^(ph) the entropy obtained from the reduced density pu for a 
given hexagon, and similarly for Sp and St- Due to the overlap of neighbouring hexagons, 
we must subtract the contribution of pairs of neighbouring triangles to avoid overcounting. 
Similarly, we must add back the contribution of single triangles to account for the overlap 
of neighbouring pairs. We assume, moreover, that the reduced density pu is the same 
for each hexagon, i.e. that the distribution is translationally invariant. The variational 
principle will then be applied to this common reduced density matrix associated with the 
hexagons, which we write as: 

Pq{°1, 0"2, 03, 04, 05, 06, Zl, Z 2 , Z 3 , #4, Z 5 , Zq) (4.7) 

with the notation of Fig. 4 for the a and z variable on a hexagon. In particular, we choose a 
numbering such that odd indices refer to triangles of subset A and even indices to triangles 
of subset B. The above reduced density matrix represents the probability for one hexagon 
to have fixed values of a and z. It is non-zero only for the allowed sets of a and z satisfying 
the two basic folding rules, and it is normalised to one. We also introduce the two-point 
and one-point functions to account for Sp and St in (f4.6|): 

p 2 ((Tl, 02, Zi, Z 2 ) = - E [Pe(o*l> °2, 03, O4, (T5, (T 6 , Z t , Z 2 , Z 3 , Z 4 , Z 5 , Z 6 ) 

< T 3' <T 4' CT 5' CT 6 

+ P6(°"37 02, 0*1) 04, ff 5, °~6i z 3i z 2, z l, &4) Z5, Zq) 

+ P6(<J3, 04, cri, 2 , 05, 06, Z3, Z4, Zi, Z 2 , Z5, Zq) (4.8) 

+ Pe(03, 04, 05, 02, 01, 06, Z%, Z4, Z 5 , Z 2 , Zi, Zq) 

+ Pe(03, 04, 05, 06, 01, 02, Z3, Z4, Z 5 , Zq, Z\, Z 2 ) 

+ Pe(01, 06, 03, 04, 05, 02, Zi, Zq, Z3, Z±, Z5, Z 2 )\ 
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and 

Pia(cti, zi) = Y P&li a 2, z u z 2 ) , 

(4.9) 

Pib{oi, z%) = Y P2(<7l, °2, z u z 2 ) . 

cri, z-y 

The two-point function p 2 is constructed to be invariant under rotations of the hexagon 
by 27r/3, i.e. by those rotations which leave the two subsets A and B invariant. Using the 
CVM equation below, this will in turn imply that the six-point function p§ itself has this 
27r/3 rotational symmetry. In ( |L9| ) we have introduced two site density matrices, pxA an d 
Pib, corresponding to the two subsets A and B. 

Noting that the numbers Njj, Np, Nta arid Ntb of hexagons, pairs and triangles of 
the subsets A and B, respectively, satisfy Np/Nn = 3, Nta/Nh = Ntb/Nh = 1, one 
finds the following approximate CVM free energy per hexagon: 

fMWi}, {zi})) = -K Y 0*1 ^2(1 + 2z 1 z 2 )p 2 (o-i, cr 2 , z u z 2 ) + A Y Ps 

1 ' 1 1 ' 1 (4.10) 

+ Y (Peinpe) ~ 3 Y (P^ ln P^) + Y (Pia^Pia) + Y (Pifllnpis), 

{cr,z} {>,z} {°": 2 } {^i*} 



which is a functional of P6({cri}, { z i}) on ^J (by implicit use of Eqs. (|4.8|) and ( |4.9|) ), to be 
minimised with respect to Pe({o~i} ^ {zi}) . Here A is a Lagrange multiplier which ensures 
the normalisation of pe({ai}, {zi}) to one. The derivative with respect to a generic element 
of p% yields the stationarity conditions 

Pe(Wi}i {zi}) = exp[-A + y Y 0- + 2ziZi+i)criO-i+i] 

i=l,6 

x [P2(o"i, oe, «i, z 6 )p 2 (cri,a 2 , z u z 2 )p 2 (a 3 , a 2 , z 3 , z 2 ) 

P2(0"3, 0"4, ^3, ^4)102(^5, 0"4, «5> ^4)102(^5, CT 6 , Z 5 , Z 6 )] 1/2 

x [pia(o"i, zi)pib{(T2, z 2 )pia(o- 3 , 2; 3 )pis(<T 4 , 24)^1,4(0-5, 2 5 )pib(o~6, ^e)]~ 1/3 , 

(4.11) 

with the convention cr 7 = &x, z 7 = z\. 

One can solve this set of equations, with the definitions of Eqs. (|4.8|) and (|4.9|) , by 
iteration. Starting from some assumed form for p 2 and iterating the above equation, 
converges to a solution of ( |4.11| ) which is, moreover, a local minimum of the approximate 
free energy Q4.10Q . At each step, the normalisation condition is imposed by adjusting the 
Lagrange multiplier A. 

In practice, we expect that the expectation values P, O and T of Eqs. ( [OP are the 
same for negative K when calculated on the subset A only or the subset B only. This 
implies pia(o~,z) = pxb{— o~, z). It is thus more convenient to use the variable cr s t, instead 
of a, and to introduce a single one point function pf , (a st , z) = pi^(cr st , z) = pis(— cr st , z) 
which is independent of the subset A or B. Eq. ( [4 . 1 1| ) , which is a set of 96 x 4 coupled 
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equations, can be reduced to a much smaller set of equations by introducing appropriate 
ratios of six-point functions and noting that some of these ratios are identical. For instance: 



p 6 (+ + -H , + + + + ++) _ pej + ++, + + + + ++) 

Pel- - + -++, + + + + ++) ~ Pe(+ + + , + + + + ++) 

(4.12) 

= P 2 (-+,++) ( pia{+,+)pib{-,+) \ ~ 3 
p 2 (+-,++) \pia(-,+)pib(+,+)J 

We can then use an appropriate parametrisation of the six-point function in terms 
of a small number of independent ratios. For these we get a much smaller closed set of 
coupled equations which can be solved numerically by recursion. In practice, we used a 
set of 9 independent ratios. 

4-2. Results of the CVM for the 2-d (arbitrary K) and 3-d (positive K) folding problems 

Before we present our results for the 3-d folding problem at negative K, let us briefly 
recall the results which have been obtained in JTJ] , where the CVM was applied to the 
2-d folding problem (i.e. to the planar 11-vertex submodel) with arbitrary K, and those 



of |T3[ for the full 3-d folding problem in the FCC lattice with positive K. 

2nd order 1st order 

Piled up Crumpled Totally flat 



K 



P2 







K, 



C2 



K 



Fig. 8: The phase diagram in the bending rigidity K of the 2-d folding prob- 
lem. 



For the 2-d folding problem, characterised in our notation by z = +1 everywhere, the 
CVM predicts two transitions, as depicted in Fig. 8. A first order crumpling transition 
at Kc2 ~ 0.1013 separates a large positive K totally flat phase with a = +1 everywhere 
from a low \K\ crumpled phase with (a) = 0. This latter phase also has (cr a t) = and 
persists with decreasing K to negative values until Kp2 ~ —0.284, where a second order 
piling transition takes place.i For K < Kp2, the system develops a non-zero value of (cr st ), 
indicating a tendency for triangles to pile on top of each other. At K = 0, the number 
of folded states for triangles is found to grow like q N with q = y/l3/3 ~ 1.202, to be 
compared with the exact value q = v / 3r(l/3) 3//2 /27r ~ 1.208. Note that this phase diagram 
is also found by another approach to the problem using a transfer matrix formalism || . 



An exact formula for Kp2 can be found in [19]. 
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We stress here that, thanks to the duality transformation, this phase diagram applies 
to the 11-vertex model describing folding on a single tetrahedron. The only modification 
is that the energy difference between the two types of allowed folds, which is K — (-K) = 
2K in the planar problem, is now —K/3 — (-K) = 2K/3. The phase diagram of the 
tetrahedral folding problem is thus exactly the same as that of the planar folding problem 
with K — > 3K. For instance the piling transition now occurs at K ~ —3 x 0.284 = —0.852. 

For the 3-d folding problem in the FCC lattice, the CVM method has been used 
in |13| for positive K only. There it was found that, as in the 2-d first order 

crumpling transition takes place, now at a larger value Kc3 ~ 0.1855 between a large K 
totally flat phase with z = a = +1 everywhere, and a crumpled state with (z) = (a) = 0. 
A remarkable difference is that the low K crumpled state now has (cr s t) non-zero. For 
instance, at K = 0, one gets (cr st ) ~ 0.87456. As mentioned at the beginning of this 
section, this indicates a strong preference for the lattice to wrap on octahedra. This non- 
zero value of the staggered a variable is not surprising in the sense that the 3-d folding 
problem in the FCC lattice may be formulated approximately as a 2-d folding problem with 
a staggered magnetic field (see flT2f ). Furthermore, a comparison, at K = 0, of the entropy 
factor q for the full 3-d folding problem (estimated to be 1.428 by the CVM and 1.43 
by a transfer matrix approach ||12|] ) with that of the single octahedron folding problem 
(q = v2 ~ 1.414) shows that the wrapping of octahedra is the entropically dominant 
folding process at zero bending rigidity. 

4-3. 3-d folding at negative K 

We now present our results for the CVM applied to the full 3-d folding problem in the 
FCC lattice, with negative K. Our main result is the existence of two separate geometrical 
transitions with increasing \K\: a first discontinuous transition at K = Kqz. ~ —0.294 
separates a phase where the triangular lattice is preferentially wrapped around octahedra 
from a phase where it is preferentially wrapped around tetrahedra. A second continuous 
transition at K = Kp3 ~ —0.852 separates this latter phase from a phase of complete 
folding of the lattice on top of a single triangle. 

The first sign of these transitions can be seen directly by looking at the evolution with 
K of the relative proportions of the four types of folds averaged over the folded states of 
the lattice. These proportions are displayed in Fig. 9. At K = 0, the fraction of acute folds 
and no- folds is identical, as a consequence of duality. Together they represent, however, 
only 5% of all folds and the statistics is dominated by complete and obtuse folds, indicating 
a phase of wrapping of octahedra. This domination of complete and obtuse folds persists 
in the range Kq3 < K < 0. At K = Kqs, all proportions show a discontinuity. For 
K < Kg3, the statistics is now dominated by complete and acute folds, indicating a phase 
of wrapping of tetrahedra. The proportion of the other types of folds becomes essentially 
zero for K ~ —0.6. The proportion p(K) of complete folds (respectively 1 — p(K) of 
acute folds) then increases (respectively decreases) smoothly with increasing \K\. This 
domination of complete folds is enhanced at K = Kp^, where a discontinuity in the slope 
of p(K) signals a continuous transition. 

A more detailed description is provided by examining the three order parameters O, T, 
and P defined in (f4.1|) , whose values, displayed in Fig. 10, clearly indicate three separate 
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Fig. 9: The average proportions of the four types of folds, as a function of K. 



phases. At low \K\ the system is preferentially wrapped around octahedra, i.e. is in a 
phase with a non-zero value of the octahedral order parameter O, while the other two 
order parameters T and P remain zero. For intermediate \K\, the system is preferentially 
wrapped around tetrahedra, i.e. is in a phase with a non-zero value of the tetrahedral 
parameter T, while the two other order parameters vanish. The geometrical transition at 
K = Kq3 ~ —0.294 between these two phases is first order. At large values of \K\, the 
three order parameters are non-zero with T ~ +1 and O ~ P. This clearly indicates a 
phase where the triangles are preferentially piled up on top of each other. In practice, the 
tetrahedral order parameter becomes indistinguishable from its saturated value T = +1 
for \K\ > 0.6, so that at large \K\, the 96-vertex model becomes indistinguishable from 
the pure tetrahedral 11-vertex model (defined by a st = z), which satisfies O = P. It 
is therefore not surprising to recover the second order piling transition of the 11-vertex 
model, which, as explained above, occurs at K PS ~ — 3 x 0.284 = —0.852. 

The final phase diagram of the 3-d folding problem of the triangular lattice in the 
FCC lattice is shown in Fig. 11. When compared with that of the 2-d planar folding 
problem, displayed in Fig. 8, we see that the crumpled phase has been replaced by two 
geometrical octahedral and tetrahedral phases, separated by a first order transition. Note 
that the octahedral phase persists at positive K until the transition point Kc3 where the 
totally flat phase is recovered. The passage from the totally flat phase to the piled up 
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Fig. 10: The octahedral order parameter O, tetrahedral order parameter T 
and planar order parameter P as a function of the bending rigidity K. 



2nd order 



1st order 
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Fig. 11: The phase diagram in the bending rigidity K of the 3-d folding 
problem in the FCC lattice. 

phase, which involved two transitions in the 2-d case, now requires three transitions. 
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It is interesting to compare our model with the Ashkin- Teller model [20] defined in 
terms of two unconstrained Z 2 variables cr st and z with pair energy: 



E = -Ko{cr st )i(a st )2 ~ K P ziz 2 - K T (za st )i{za st )2 



(4.13) 



involving three energy parameters Ko,p,t- Note that, despite the apparent asymmetry, 
the three Z 2 variables cr st , z and za s t are actually exactly on the same footing: for any two 
of them, the third one is simply the product of these two. Our energy ( |2.5| ) corresponds 
to a trivial case (Ko, Kp, Kt) = (—K/3,0, —2K/3) where one of the coupling constants 
vanishes. The unconstrained system then translates into two independent Ising models 
with two continuous transitions with increasing \K\ describing the successive symmetry 
breakings of first (o~ s tz) at if = —3K c /2 and then of (cr s t) at K = —3K C , while (z) = 
(a s t) (za s t) for all K. Here K c = —(1/2) ln(2 — v^3) — 0.658 is the critical Ising coupling on 
the honeycomb lattice. In particular, at low \K\, all order parameters vanish. We see that 
the constraints (|2.1|) and ( |2.3|) on the cr st and z variables modify this picture by replacing 
the first continuous transition by a discontinuous one and by giving rise to a low \K\ phase 
where (cr st ) is non-zero. In this phase, the origin of the symmetry breaking in the variable 
cr st is clearly entropic. Indeed it allows for a gain in entropy from the z variable as the 
second basic folding rule is less constraining when the a variables alternate on the lattice. 
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Fig. 12: The two-point functions 02, T2 and P2 and the energy e, free energy 
/ and entropy s = e — f per hexagon, as a function of the negative bending 
rigidity K. 



Since we can extract the two-point function from the CVM, we have access to the 
nearest- neighbours correlations 02, T2 and P2 defined in ( |4.2|) , as well as to the average 
energy per hexagon: 

e = K (02 + 2T2). (4.14) 
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We also have access to the free energy / per hexagon and to the entropy part s = e — f. 
These quantities are displayed in Fig. 12 as a function of K. 

4-4- The transfer matrix approach 

The CVM gives only an approximation to the exact statistics of the problem and it is 
difficult to estimate quantitatively possible systematic errors of the method. In particular, 
due to the truncation of the free energy functional to finite clusters, we cannot even ensure 
that the CVM free energy is bounded from below by the true free energy. It is thus 
important to compare the CVM results with those of other methods. In the case of 2- 
d folding, a formulation of the problem was given in terms which allows for the use of 
a transfer matrix method 0, leading to a very satisfactory comparison with the CVM 
approach. The same method was generalised for the 3-d folding problem at K = in |E| . 
In this section, we present the results of a transfer matrix study of the 3-d folding problem 
at negative bending rigidity K. 

In the transfer matrix formalism, the folding partition function Zl(K) of an infinite 
strip with a finite width L and with rigidity K is given by 

Z L {K) = lim [Tr(Ti(.ftT) n )] 1//n , (4.15) 

n — >oo 

where Tl(K) is the row to row transfer matrix of the system. The state of a row is specified 
by the 2 x 2L triangle variables Oi and z% in the row. We use free boundary conditions 
which allow folds to escape from the strip at its boundary. In the folding problem, a 
transfer matrix element from the lower row {a, z} to the upper row {a', z'} 



T L {a,z},{a',z'} 



// 


z / 
' V ' 






z ' / 

I 2L-2 / 
a 2L-2// 
X Z 2L-1 




z 2 / 

° 2 / Z 

/ 3 

/°3 




' L-l * 


z 2L-2/ 
°2L-2 / 

/ z 2L-1 



(4.16) 



is non-zero ijf the two basic folding rules (|2.1] ) and ( |2.3|) are satisfied around each of the 
vertices j = 1, . . . , L — 1; if so, it is equal to exp(— E), where E is the folding energy ( |2.5|) , 
proportional to K, for all pairs of neighbouring triangles in the two rows. This method 
has been extensively described in [12| for the 3-d folding without rigidity (K = 0). The 
question of the rigidity has been addressed in for the 2-d folding problem. The reader 
is referred to these articles for further details on the method. 

As usual, we must extract the largest eigenvalue (and its eigenvector) of the matrix 
Tl(K). The main difficulty is that the size of this matrix increases exponentially with L 
like 4 2L x 4 2L . Also, due to the chirality in the representation of the triangular lattice as a 
square lattice in ( f!.16|) , the transfer matrix is not directly symmetric, which requires the 
use of both left and right eigenvectors. Fortunately, the vanishing of all elements which 
do not satisfy the basic folding rules makes the matrix very sparse. With the Cray C94 
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of CEA-Grenoble, we have been able to reach L = 6 with one Gigabyte of memory. By 
comparison, in 2-d folding, the size of the matrix was 2 2L x 2 2L and the width L could be 
doubled. These computations are exact, in the sense that the precision is limited by the 
numerical precision of the computer; no approximations are made with this method. Of 
course, the infinite lattice folding is recovered only in the limit where L — > oo. 

Strictly speaking, an infinite strip with finite width is a one-dimensional system which 
has no phase transition and all order parameters vanish due to the global symmetry under 
reversal of all spins. The order parameters, however, can be estimated by weakly breaking 
the symmetry by first fixing the (a, z) variables of, say, the most central triangle of a row 
i = L and by averaging each order parameter on the other triangles. In practice, it is useful 
to restrict oneself to triangles j with \j — i\ > 1, i.e. to triangles which are not nearest 
neighbours of i. In fact, the contribution of neighbouring triangles is less significant because 
the short range correlations are always large even in a disordered phase. To summarise, 
we measure Xl = 1/(21/ — 3) J2j-\j-L\>i( XLX j) as an estimate for the order parameter 
X = (x) to be considered (x = a st , za st or z, i.e. X = O, T or P). 

Using the Pade-Shanks transformation (see Ref. f2~l]1 ), the finite size effects can be 
reduced. Given three values X4, X 5 , X 6 obtained for sizes L = 4,5,6, and supposing the 
convergence is exactly exponential, the extrapolated value is X^ = (X^Xq — X§)/(X± + 
Xq — 2X5). 




Fig. 13: The results of the transfer matrix method for the tetrahedral order 
parameter T as a function of (negative) K. The figure displays the results for 
strips of size L = 3, 4, 5, 6 as well as the extrapolation to infinite L (see text). 
By comparison, the corresponding result of the CVM is also indicated. 

In Fig. 13, the tetrahedral order parameter T, estimated as described above, is plotted 
as a function of the negative rigidity K. For finite width L = 6, the first order transition 
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predicted by the CVM is smoothed by finite size effects. The extrapolation L — > oo, 
however, clearly shows a rapid change around K = —0.3. Unfortunately, the maximum 
size L = 6 is too small to allow a precise determination of Kq% or even the order of the 
transition. We can still estimate Kq3 roughly by considering the point where all curves at 
finite L seem to intersect. This leads to Kqz ~ —0.3, in good agreement with the CVM 
estimate K G3 ~ -0.294. 

For the other order parameters (P and O), the curves are smoother and the predicted 
transitions less visible. Still, the variations of these order parameters with K are in agree- 
ment with the infinite system limit expected from the CVM results. For instance, we see 
an intermediate \K\ regime with a low value of O surrounded by a small and large \K\ 
regime with a significantly larger value of O, as expected in the presence of an intermediate 
tetrahedral phase. 

Clearly, the exact transfer matrix method is more difficult in the 3-d folding case than 
in the 2-d folding case, due to the higher number of states. In the 2-d case, the first signal 
of the continuous piling transition at Kp 2 did not occur until system width L = 8. Such a 
large system size is impractical in the 3-d problem. In this context we can conclude that 
the transfer matrix results do not contradict in any way those of the CVM, which in turn 
appears to be a much more powerful approach. 

5. External fields 

In this section, we briefly discuss the consequence of adding external fields coupled to 
our Z<i spin variables. The results of this section are again obtained in the framework of 
the CVM. We add to the pair energy ( |2.5| ) the following external field contribution: 

E H = -H Ost - H T za st - H P z , (5.1) 

with three external fields Ho, Ht and Hp coupled to the order parameters O = (cr st ), 
T = (za st ) and P = (z), respectively. 

For simplicity, we will discuss in the following the consequence on the negative K 
phase diagram of having one of these fields non-zero, while the two others are kept equal 
to zero. 

As explained in Section 4.3, in the absence of external field, the system was found 
to have three phases for negative K. The nature of each phase can be characterised by 
specifying which of the three order parameters O, T and P are zero or not. 

(i) The low \K\ octahedral phase has O ^ 0, P = T = 0. 

(ii) The intermediate \K\ tetrahedral phase has O = 0, P = and T ^ 0. 

(iii) The large \K\ piled up phase has O ^ 0, P ^ and T ^ 0. 

We find that the nature of a phase is unchanged by the addition of a (small enough) 
external field Hx if the order parameter X was non-zero in this phase. Here X stands for 
O, T or P. On the other hand, if X happened to be zero in a phase, then the action of 
Hx is to make all order parameters non-zero. 

For instance, the piled up phase, which has its three order parameters non-zero, per- 
sists in the presence of external fields. In the presence of a non-zero value of Ho, the octa- 
hedral phase is preserved but the continuous piling transition is suppressed. The original 
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tetrahedral and piled up phase merge into a single phase with non-zero order parameters 
O, T and P. This phase can be viewed as a tetrahedral phase with a smoothly increas- 
ing piled up order with increasing \K\. The first order transition, now directly from the 
octahedral phase to this phase, is slightly shifted to larger values of \K\. In the presence 
of a non-zero value of Ht, the tetrahedral and piled up phases are not modified and the 
second order transition is preserved at its original value of K. This is not surprising since 
T was already saturated to T ~ 1 in this regime. On the other hand, the system develops 
non-zero order parameters T and P in the low \K\ region, indicating the appearance of a 
piled up order in the original octahedral phase. The first order transition from this phase 
to the tetrahedral phase is preserved (for small enough Ht at least), although slightly 
shifted to smaller values of \K\. Finally, in the presence of a non-zero value of Hp, the 
three order parameters become non-zero over the whole negative K range. The second 
order transition is suppressed while the first order transition now takes place between two 
phases of the same nature, discriminated only by which fraction of acute or obtuse folds 
is larger. 

The above discussion is illustrated in Fig. 14 where the three order parameters O, P 
and T and represented as a function of K in the presence of an applied external field Ho, 
Ht or Hp. 



6. Discussion 

In this paper, we have studied the folding of the 2-dimensional regular triangular 
lattice in the 3-dimensional regular Face-Centred Cubic lattice, a problem which can be 
expressed as a 96- vertex model on the triangular lattice. We considered the case of negative 
bending rigidity K, where folds are favoured. A special role is played by three classes of 
elementary objects in the FCC, namely elementary octahedra, tetrahedra and triangles. 
Indeed, instead of randomly wandering over all sites of the FCC lattice, the system prefers 
to stay trapped on one class of these objects: 

- At low \K\, the system is preferentially wrapped on octahedra, characterised by obtuse 
109° degree folds. 

- At intermediate \K\, the system is preferentially wrapped on tetrahedra, with acute 
71° folds. 

- At large \K\, the system is preferentially piled onto single triangles, making complete 
0° folds. 

This leads to three distinct phases, characterised by the symmetry breakings of appro- 
priate order parameters. These phases are separated by two transitions, a discontinuous 
one between the octahedral and the tetrahedral phase and a continuous one between the 
tetrahedral and the piled up phase. This study completes the previous one at positive 
K which found a first order crumpling transition with a large positive K completely flat 
phase. 

Our results were obtained in the framework of the CVM approximation and further 
corroborated by a transfer matrix approach. At this stage, one may wonder whether the 
particular form Q2.4] ) of the bending energy is crucial for our results. We expect that, as 
long as the energy is a monotonic function of the angle of the fold which, moreover, favours 
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Fig. 14: The variation with the negative bending rigidity K of the three order 
parameters O, T and P in the presence of a small applied external field. The 
first (respectively second and third) column corresponds to the application 
of a positive external field Ho (respectively Ht and Hp). The dashed verti- 
cal line indicates the location of the first order geometrical transition in the 
absence of external fields. 

folding, the above scenario will continue to hold. One may also question the validity of the 
CVM approximation in a very constrained problem where folds are, moreover, expected 
to propagate throughout the entire lattice, far beyond the size of the basic cluster of 
the CVM. Actually, as explained in || for the 2-d folding problem, the necessity for 
folds to propagate through the lattice is entirely encoded in the particular form of the 
96 local hexagonal configurations and the matching condition for the folds on links which 
are common to two overlapping hexagons. These two ingredients are properly taken into 
account by the CVM. We thus believe that this method does indeed account for the 
propagation of folds and therefore gives reliable results concerning the phase diagram and 
the transitions separating distinct phases. Nevertheless, as for any mean field approach, 
the method gives only trivial critical exponents for the continuous piling transition, whose 
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true universality class is therefore still to be determined. A CVM approach with larger 
clusters would not help to resolve this issue and we believe that the restriction to hexagonal 
clusters is probably the most efficient choice in terms of simplicity. On the other hand, a 
transfer matrix method, which is exact for strips of finite size L does not allow us to reach 
large enough system sizes to allow a reliable continuum (L — > oo) extrapolation. It is thus 
a challenge to find alternative methods for a more refined study of this problem. 

Finally, the particular sequence of folding involving octahedra and tetrahedra is clearly 
intimately linked to the structure of the FCC target space where the 3-d rotational sym- 
metry has been explicitly broken to a smaller (cubic) symmetry. On the other hand, the 
piled up phase with alternating normal vectors on neighbouring triangles is likely to exist 
even if the target space is continuous. One may wonder whether the folding of the triangu- 
lar lattice in a continuous 3-d space also gives rise to intermediate phases with particular 
broken symmetries, not necessarily related of course to the FCC cubic symmetry. 
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